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We use a combination of perturbation theory and numerical techniques to study the equihbration 
of two interacting fields which are initially at thermal equilibrium at different temperatures. Using 
standard rules of quantum field theory, we examine the master equations that describe the time 
evolution of the distribution functions for the two coupled systems. By making a few reasonable 
assumptions we reduce the resulting coupled integral/differential equations to a pair of differential 
equations that can be solved numerically. Our results show with good accuracy how these coupled 
I ■ systems approach a common equilibrium temperature. 
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, Non-equilibrium physics is a topic of increasing interest in physics today. Such effects are expected to play an 

' important role in many contexts, such as the formation of the quark gluon plasma |l[ and electroweak baryogenesis 

d [H . Compared to the equiUbrium case, the details of non-equilibrium processes are difficult to analyze, and a variety 

5! of different approaches have been used: analysis of the equations of motion |^ , renormalization group techniques [|| , 
kinetic and transport |6| theory, and thermo field dynamics Q. In this article we develop a new approach, and 
apply it to a simple interacting model of two fields initially each in thermal equilibrium but at different temperatures. 

I , We find we can reliably track the path towards equilibration of this system using perturbation theory applied to the 

^ master equation describing the distribution functions of these fields. The analysis allows us to obtain quantitative 

00 information about this non-equilibrium process directly from the microscopic field theoretic description of the physical 

0^ system. The method can readily be generalized to more complicated, and physically relevant, systems. 

, The partition function describing the system we are considering is given by, 

2 : Z = j d0exp{* j dx[]^{d^(^af - + \{d^(^bf - m^l ~ ^-\<PaM^]} (1) 

, where a and b are the two field types. We will use the Keldysh representation of the closed time path formulation of 
finite temperature field theory Q, so that indices {1, 2} indicate the top/bottom branch of the time contour. 
^ , The Schwinger-Dyson equation for the two-point function for field a has the form , 

[{d-dy,-{d-d)y\SU^,y)= Jdz [S'^,ix,z)^Uz,y)-SU^,z)^^,^iz,y) 

-S?i(x,z)5«2(z,2/) + S?2(x,z)5«2(z,2/)] (2) 

K> . . . 

: and similarly for system b. Writing all functions in terms of an average X = ^{x+y) and a relative s = x—y ; s' = x — z 
d ' coordinate, and dropping terms (s — s') and s' relative to X, taking the transform J dse^^'^^ of (||) yields 

- 2ik^^^S^,{X, k) = S^^iX, k)T.l^{X, k) - k)SUX, k) (3) 

We now make the ansatz for the propagators, 

ZTT 

SUX, k) = -[5{k^ - u;t)n''iko, X) + <5(fco + uol){l + n'^i-ko, X))] 
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S^,{X, k) = -[5{k^ - Lol){l + n-{k,,X)) + 5(fco + c.^)n'^(-fco, X)] 

^k 

where LOa{k) = Assuming now spatial homogeneity, equating the positive frequency parts in @ gives 

dna{u}^,t) 



dt 

where n"(a;^,X) = and 



= -n,K,<)rrfK) + [1 + n,K,t)]r,K) (4) 



^^("'^") = 2if^"^' = 2^^^! (5) 



2 



are, respectively, rates associated with processes which increase/decrease the number of particles present of that 
species. These rates are complicated functionals of the distribution functions, and so (^) is very difficult to solve in 
general. However, matters become simpler in the case that each distribution function is, at a given time, close to a 
thermal equilibrium distribution, as we will now describe. 

We assume that the two fields start out in (different) thermal equilibrium states with initial conditions 



n^(fco,0) 



1 



1 ' 



— a.b 



(6) 



We further assume that the interactions are described by the two-loop sunset self-energy. For given thermal equilibrium 
distribution functions, a detailed calculation reveals that: 



167r2 



dp—: 



dx 



[(1 -f n{u;Pa))C\ko, k,iu;,p) + n{iu;f3a)C\ko, k, -u;;,p)] 



(7) 



S?2(fco,A:) 



„2 /•! 



dp ■ 



dx 



J-i 



[n(o.;/3,)/:''(c.«,p, fco, ^) + (1 + n(c.;/3,))/:^(-c.;,p, fco, fc)] 



(8) 



where uJx{k) = ^yk^ + rn^ and the function £^{K,P) is given by, 

52[l + n^'{(fc0-/)/34] 



C\K,P) 



8'K\k-p\Pb 

+ OiiK^Pf-AniDln 



X <j e{-{K-Pf)\n 
sinh(r+/3b/2) 



1 - exp(-r+/3b) 



1 - exp(-r_/3;,) 



sinh(r_/3b/2) 



(9) 



where 



r± EE ^ ( |fc-plVl-4m2/(i^_p)2±(fcO_pO- 



|fc-pl yjk'^+p^ - 2kpx, P^^pI-p^. 
It is important to note that for the initial conditions (||), the KMS condition 



(10) 



(11) 



is satisfied if and only if Pa = Pb] i-c the combined system is in thermal equilibrium. 

We take the following approach to solving (^. We begin with initial conditions (^) which describe the two fields 
at two different (inverse) temperatures j3a and /?(,. We assume that the equilibration time for each separate system is 
very small compared to the time scale of equilibration for the coupled system, so that the distribution function for 
each system remains close to a thermal one as the systems evolve. Now, near t — Q, the distribution functions which 
appear within the rates in (^ can be replaced by their thermal values, and (^ has the solution 



+ C(u;^)exp(-[rd-r,]t) 



(12) 



with C{u)1) determined from the initial condition. Note that if [3a = /Sh, the KMS condition (|ll| ) is satisfied, and 
substituting into (^ we find C(w^) = 0, so that na{oj'^,t) = na(a;^,0); this simply means that the combined system 
stays at equilibrium, as expected. Thus, for a suitably short time interval 5t <^ l/p^ — Fi], ( |l2|) gives the approximate 
time evolution of each distribution function. To integrate (^ to finite time intervals, we first let the systems evolve 
according to (12) for a time 5t, and the subsystems are assumed to come separately to thermal equilibrimn. We then 
fit the distribution functions resulting from (|l2[) to equilibrium distribution functions of the form: 



nx{ko, \k\) 



(13) 
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by adjusting both the parameters characterizing the inverse temperature (3x and a (temperature dependent) mass 
rux for each subsystem. The thermal distributions resulting from this best-fit procedure are then used as the initial 
conditions for the next time iteration. This results in a second evolved set of distribution functions which again are fit 
to thermal distributions. This process is continued until (hopefully) the best fit temperatures and masses, and hence 
the new distribution functions, change little from one time step to the next. 

We present some results of this procedure for initial temperatures of 0.7 and 0.75 and initial masses of 0.6 and 0.65 
for, respectively, the a and h fields. The resulting evolution, for time steps of 50 units, for the inverse temperature 
and masses appears Fig. dl^) - note that in both cases a common equilibrium is reached. 




4000 



FIG. 1: Evolution of the inverse temperature (left) and the mass (right). The crossed points represent the a field, and the 
boxed points represent the b field. 



Using these distributions we can follow, for each species, the evolution of the particle number N{t) and energy E(t): 

N{t) = / k^dk n{k,f3,m) 



1 foo 

E{t) = —FT k^ dk n{k,l3,m) (14) 

Jo 

Again for time steps of 50 units, these results appear in Fig. (^) - to a good approximation the total particle number 
and the total energy remains constant, although small numerical deviations are observed. 



Particle number evolution 




FIG. 2: Evolution of the particle number (left) and the energy (right). The crossed points represent the a field, the boxed 
points represent the b field, and the triangulated points of the top line represent the sum. 

The preceding results were obtained using a basic time step of 50 units. We can check the consistency of the results 
by following the evolution of the system for different basic time steps and then comparing the results at common 
times. Results are shown in the following table for the evolution of (3 for field a using basic time steps of 50, 100, 200, 
400, and 800 units. 
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TABLE I: Comparison of the evolution of the inverse temperature of a for basic time steps of 50, 100, 200, 400, and 800 units. 

We see that the results amongst the different time steps agree to a good approximation at the corresponding times. 
However, as is expected, for larger time steps some disagreement starts to emerge. 

In conclusion, we have shown that equilibrium is achieved in the simple model we have considered. Starting from 
two separate systems at temperatures that are close enough together so that the combined system is never very 
far from equilibrium, we fit to a series of equilibrium distributions using time steps that are small compared with 
the equilibration time. The result is that the combined system reaches equilibrium at a temperature between the 
temperatures of the two original systems. In addition, we have shown that both particle number and energy are 
conserved. It would be interesting to try to relax the assumption that the initial temperatures of the two initial 
systems are close together. Allowing for a larger difference in initial temperatures would create a more highly non- 
equilibrium state in the early stages of the evolution. Such a situation might be handled by using a more complex 
ansatz for the distribution functions of the evolving system. 
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